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Abstract 



O 

c/3 It is shown that the dispersion and attenuation functions in a hnear viscoelas- 

'^ tic medium with a positive relaxation spectrum can be expressed in terms of 
"S a positive measure. Both functions have a subhnear growth rate at very high 
. frequencies. In the case of power law attenuation positive relaxation spee- 
ch trum ensures finite propagation speed. For more general attenuation func- 
S tions the requirement of finite propagation speed imposes a more stringent 
i-rt condition on the high-frequency behavior of attenuation. It is demonstrated 
fl that superlinear power law frequency dependence of attenuation is incom- 
patible with finite speed of propagation and with the assumption of positive 
relaxation spectrum. 
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Laplace transform; 
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Laplace transform; 
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homogeneous distribution; 
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Caputo fractional derivative 
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< lim^^a f/g < 00; 
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1. Introduction. 

Viscoelastic media of the Newtonian, hereditary or combined type are 
known to have positive relaxation spectra [U EJ [3l HI |5]. Any network of 
springs and dashpots combined in parallel and serially has this property [6]. 
This is also true of many theoretical models of polymers (the Rouse model 
[7], Zimm's reptation-based models [S], Generalized Gaussian models [H]), 
rubber [TUl [H] , metals [12] and rocks [T3] . Universality of this property has 
led to same attempts at explaining it by a general thermodynamic or "fading 
memory" principle |21[Illll])- The same property is shared by some (but not 
all) dielectric media [15]. 

Positive relaxation spectrum has some consequences for wave dispersion 
and attenuation in such media. We shall use the term "admissible" to refer 
to the attenuation and dispersion functions which are compatible with posi- 
tive relaxation time spectrum. It will be shown that every admissible pair of 
attenuation and dispersion functions has an integral representation in terms 
of a positive Radon measure, which we call the dispersion-attenuation mea- 
sure. The measure is arbitrary except for a very mild growth condition. The 
two expressions taken together can be considered as a dispersion relation in 
parametric form. 

In contrast the acoustic Kramers- Kronig (K-K) dispersion relations [16] 
are non-local. The K-K dispersion relations express the dispersion function 
in terms of the attenuation function or conversely. This presupposes that 
one of these functions (usually the dispersion function) is accurately known 



and and is admissible. On the other hand, substituting any positive measure 
in the parametric dispersion relation yields an admissible pair of dispersion 
and attenuation functions. 

There is no a priori argument in favor of the K-K dispersion relations 
in acoustics. In electromagnetic theory the K-K relations follow from the 
causality of the time-domain kernel appearing in the constitutive equations 
of a dispersive dielectric medium. In acoustics the K-K relations are derived 
from an ad hoc assumption about the analytic properties of the wave number. 
It is namely assumed that the wave number is the Fourier transform of a 
causal function (or distribution) L. This assumption is not justified unless 
a physical meaning is given to the causal function L. Attempts to give a 
physical meaning to the function L in terms of an ad hoc wave equation are 
discussed in Sec. [6j The wave equation is however not consistent with the 
constitutive equations of viscoelasticity unless an approximation are made. 

It will be also shown that in media with positive relaxation spectrum 
the frequency dependence of the attenuation function Aiiu) in the high fre- 
quency range is sublinear: Aiiu) = o[u]. This fact raises some problems 
which we also try to address in this paper. In the case of power law attenua- 
tion the attenuation and dispersion are proportional to a power of frequency 
\u}\°'. Numerous experiments in acoustics indicate that the power law ac- 
curately represents the frequency dependence of dispersion and attenuation 
over several decades of frequency. If the power law extends to arbitrarily 
large frequencies then viscoelastic theory based on positive relaxation spec- 
trum implies that the exponent of the power law satisfies the inequalities 
< a < 1. The same inequality follows from the requirement of finite speed 
of propagation. 

Investigations of creep and relaxation in viscoelastic materials always sup- 
port the assumption of positive relaxation spectrum (e.g. [I7l[T2] for creep in 
metals, |T3] for the upper mantle) and therefore models derived from constitu- 
tive relations exhibit sublinear attenuation and dispersion at high frequencies. 
In particular this applies to the Cole-Cole, Havriliak-Negami, Cole-Davidson 
and Kohlrausch- Williams- Watts relaxation laws commonly applied in phe- 
nomenological rock mechanics, polymer rheology p!8], bio-tissue mechanics 
(e.g. for bone collagen [19]) as well as for ionic glasses [20j. It also applies 
to the most common constitutive model of bio-tissue based on power law 
creep [21] or to the soft-glassy model of cytoskeleton mechanics [221 123] . 

Sublinear power law attenuation (a < 1) has also been reported in exper- 
imental investigations involving PVC and lava samples pi]. Experimental 



ultrasound investigations of numerous materials point however to higher val- 
ues of the exponent in the power law attenuation. For example in ultrasound 
investigations of soft tissues the exponent varies between 1 and 1.5, while 
in some viscoelastic fluids such as castor oil it lies between 1.5 and 2. In 
multi-walled carbon nanotubes the exponent lies above 1.1 [25J. We shall 
refer to this case as superlinear frequency dependence. Typical values of the 
power law exponent in medical applications using ultrasound transducers are 
a = 1.3 in bovine liver for 1 -=- 100 MHz, a = 1 -^ 2 in human myocardium 
and other bio-tissues [26]. Values in the range 1 -i- 2 are observed at lower 
frequencies in aromatic polyurethanes [27]. Nearly linear frequency depen- 
dence of attenuation is well documented in seismology [2B] • Approximately 
linear frequency dependence of attenuation has been observed in geological 
materials in the range 140 Hz to 2.2 MHz. 

Several papers have been devoted to a theoretical justification of the 
superlinear dispersion-attenuation models [291 EOl EH [321 |33l JMl |33]. In 
order to resolve some problems Chen and Holm [32] suggested to add a 
fractional Laplacian of the velocity field of order y, where < y < 2, to the 
usual Laplacian of the displacement field in the equations of motion. Their 
approach still leads to an unbounded sound speed for superlinear attenuation 
and adds a new problem: their wave equation cannot be derived from a 
viscoelastic constitutive equation. The underlying physical model of Chen 
and Holm's equation is for the time being unclear. 

The main factor affecting viscoelastic behavior and viscoelastic wave 
propagation is relaxation represented by the relaxation modulus. This re- 
mains true for fractal media such as polymers and bio-tissues. In contrast 
to anomalous diffusion, viscoelasticity is incompatible with fractional Lapla- 
cians. 

Superlinear power law attenuation is incompatible with a finite upper 
limit on propagation speed. This also applies to more general attenuation 
functions with a linear or superlinear asymptotic growth in the high frequency 
range. In fact we shall see that the boundary between finite and infinite 
speed of propagation corresponds to asymptotic attenuation ~ aco'/ln(c<;). 
An abnormal feature of wave propagation in media with superlinear power 
laws (and more generally in media with superlinear asymptotic frequency 
dependence) is appearance of precursors. The precursors extend to infinity 
and thus the speed of propagation of disturbances is infinite. 

It is a challenging problem how to explain the incompatibility between 
the theory of viscoelasticity (or finite propagation speed) and experiment in 



the case of superlinear attenuation. It seems likely that the attenuation ob- 
served at ultrasound frequencies should be associated with an intermediate 
frequency range in which the attenuation function behaves in a significantly 
different way from its asymptotic behavior. However even in this case prob- 
lems remain. A superlinear power law cannot be strictly valid in any finite 
range of frequencies. Since the local exponent ai{u) := ln^i(a;)/lnci; is an- 
alytic for w > 1, it cannot be constant over a finite frequency range without 
being constant over the entire real line. Hence a local superlinear power law 
entails superlinear asymptotic behavior which is incompatible with the re- 
quirements of positive relaxation time spectrum and finite propagation speed. 

2. Constitutive assumptions and basic definitions. 

In viscoelasticity the relaxation modulus G, defined by the constitutive 
stress-strain relation 

a{t) = f G{t- s) e{s) ds + ( e{t) (1) 

Jo 

is assumed to have positive relaxation spectral density h. The latter state- 
ment means that for every t > 



G{t) 



POO 

Go+ / e-*"/i(r)dln(r), t>0 (2) 

Jo 

where r = 1/r is the inverse of the relaxation time and h(r) > 0. Eq. pi) rep- 
resents a superposition of a continuum of Debye elements. For mathematical 
convenience eq. (|2| will be replaced by a more general equation 



G{t) = / e-*Xdr), t>0 (3) 

J[0,oo[ 

where /x is a positive measure: fi{[a,b]) > for every interval [a,b] of the 
positive real axis. As indicated in the subscript of the integral sign the range 
of integration is the set of reals satisfying the inequality < r < oo. In 
general the measure /i({0}) of the one-point set {0} is finite and equal to the 
equilibrium modulus Goo '■= limi_>.oo G'(t). An additional assumption 



fi{dr) 

[0,oo[ 1 + ^ 



< cx) (4) 



ensures that the function G is integrable over the interval [0,1] [S]. The 
relaxation modulus assumes a finite value (7(0) = M at if the measure /i 
has a finite mass M = /i([0, cxd[). 

In contrast to eq. ^ the integral representation ^ includes as special 
cases finite spectra of relaxation times corresponding to superpositions of a 
finite number of Debye elements 



N 



G(t) = ^c„e-'-"*, c„>0,r„>0 forn = l,...,iV (5) 

n=l 

(Prony sums), infinite discrete spectra corresponding to Dirichlet series as 
well as discrete spectra embedded in continuous spectra. However the main 
advantage of (tsl) over ^ is the availability of a very powerful mathematical 
theory which ensures logical equivalence of certain statements about material 
response functions. 

In particular a function satisfying ([s]) is completely monotone. A function 
G{t) is said to be completely monotone if it continuously differentiable to 
every order and 

(— l)"'D"G'(t) > for all non-negative integers n and t > (6) 

Bernstein's theorem [HSl [H7] asserts that eq. ([S]) is equivalent to (|6|. There 
is no such simple characterization of functions which have the integral rep- 
resentation ([2]). Hence the requirement of positive relaxation time spectrum 
(Is]) is equivalent to the complete monotonicity of the relaxation modulus 

[381 [H [391 E]. 

For us the main benefit from using (|3| instead of (|2| is the equivalence 
of ^ with a property of the dispersion and attenuation that will be ex- 
plained below. That is, certain statements about attenuation and dispersion 
functions follow from ^ or, conversely, imply that ^ does not hold. 

We now recall the definition of a Bernstein function [ID]. A differentiable 
function / on ]0, oo[ is a Bernstein function if / > and its derivative /' is 
completely monotone. A Bernstein function is non-negative, continuous on 
]0, oo[ and non-decreasing, hence it has a finite value at 0. A function / on 
[0, oo[ which has the form f{x) = x^ g{x) for some Bernstein function g is 
called a complete Bernstein function (CBF) [37] . It can be proved that every 
complete Bernstein function is a Bernstein function |37j . 

The following facts about complete Bernstein functions will be needed 
here. 
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Theorem 2.1 

A real function f on ]0, oo[ is a complete Bernstein function if (i) f has an 
analytic continuation f{z) to the complex plane cut along the negative real 
axis; (ii) /(O) > 0, (Hi) f {z) = f{z), (iv) lm.f{z) > in the upper half plane 
Iraz > 0. 



Theorem 2.1 has the foUowing coroUaries|41j: (1) If / is a CBF then /" is a 

CBF if < a < 1. 

(2) If / ^ then the function x/f{x) is a CBF. 

Every complete Bernstein function / has the integral representation: 

f{x) = a + bx + xf ^ (7) 

where a,b > and i/ is a positive measure satisfying the inequality 

i^(dr) ,„, 

TT7 < °° (8) 

]0,oo[ i +^ 

[37] . Note that the integration domain does not include the point 0. The 
contribution of would have the form z/({0}) x, hence it could be included in 
the second term on the right-hand side of (JTl). If the point is excluded then 
the constants a, b and the measure u are uniquely defined by the function /. 
If G is a completely monotone function integrable over [0, 1] and satisfying 
eq. ([3]) then 

G(p) = f ^^i^ = M{0}) + f /^(dr) 



'[o,oo[ P + r P J]o,oo[ P + r 

It follows that 

Q{p):=pGip) = Km+P [ ^ (9) 

i]o,oo[ P + r 

is a complete Bernstein function [H]. 

A function / is called a Stieltjes function if it has an integral representa- 
tion 

/W=a+/ ^ (10) 

where a > and i/ is a positive measure satisfying (|8 



3. Dispersion and attenuation. 

The Green's function Q{t, x) is defined as the solution of the problem 

pu^tt = G{t) * u^txx + S{x) 6{t) (11) 

with zero initial data. In a three-dimensional space 

G{t,x) = 



gifer 



"* T^TT dp / TT 7z^^ k dk 



where 



(271)3 ri_i^+, Q(p) -J_^ k^ + B{py 

-^ / -^e^*-^(^)^'dp (12) 

It will be assumed that G is completely monotone and that the equilibrium 
modulus Goo := linit^ooG(t)0 is positive. The first assumption entails that 
Q is a complete Bernstein function. From the corollaries of Theorem 2A we 
know that the functions Q^/^ andp/Q(p)^/^ are complete Bernstein functions. 
Hence 5 is a complete Bernstein function. The second assumption entails 

that Q(0) = limp^o pG{p) = G^o > 0, hence B{0) = 0. In view of Q 

the function B has a representation B{p) = cp + b{p) where c > 0, the 
dispersion-attenuation function b has the integral representation 

Hp):=pf ^ (14) 

J]o,oc[ P + r 

and i^ is a positive measure satisfying the inequality (|8|. A function b hav- 



ing an integral representation (14) will be called an admissible dispersion- 
attenuation function. The measure z/ is the spectral measure of the dispersion- 
attenuation function b. A dispersion-attenuation function 6 of a viscoelastic 
material with a positive relaxation spectrum is admissible. The converse 
statement is not true, as will be shown later. 

We shall now prove that limp_5.oo B{p)/p = c. The last conclusion follows 
from the fact that for p > 1 

1 1 

< 

p + r 1 + r 



and the right-hand side is integrable with respect to the measure v in view 
of eq. (|8]) . For p — )■ oo the integrand of 

z/(dr) 



]o,oo[ V^r 

tends to zero, hence, by the Lebesgue Dominated Convergence Theorem, the 
integral tends to zero as welL It follows additionally that 6(p) = o\p\ for 
p — )■ oo. We now note that limp_>oo Q(p) = Co := limi^oG'(t) is the elastic 
(instantaneous response) modulus. Hence c = I/cq, where cq := yGoJp- 



Hence eq. ( 12 ) can be recast in a more explicit form 

"ioo+e 



''^' ' 8nHrJ_,^^, Qip) ^ ^ ^ 

and the dispersion-attenuation function b{p) has sublinear growth. 

We also note that the attenuation function is non-negative in the right- 
half of the complex p-plane: 

A(p) := Reb(p) = [ IpI +^^^P ;,(dr) > (16) 

ilo,oof b + r|2 



for Rep > 0. 

The derivative 



b'ip) = M_ f p^. < M (17) 

P J]o,oc\ [p + ry p 



Moreover 

J]0,oo[ [P + "J 



On account of (tsl) the function r/(p + rY < 1/(1 + r) (for p > 1) is integrable 
with respect to the measure u and tends to zero for p — )■ oo. Hence b'{p) = o[l] 
for p —¥ CO. 

The attenuation function A and the dispersion function V{p) := — Im b{p) 
satisfy linear dispersion equations in parametric form: 

, / N f IpP -|- r Rep , , , , , 

Aip) = / ^V— nr^Kdr) (18) 

^lO,oof \P ~r i\ 



|p _|_ ^|2 



V{p) = -lmp j--^^^u{dr) (19) 

J]0,oo[ 
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The measure z/ represents the dispersion-attenuation spectrum. An ele- 
mentary dispersion-attenuation is represented by the function p/{p + r) for 
a fixed value of r. Given a CBF function 6, the ratio g{p) := b{p)/p is a 
Stieltjes function. If the dispersion-attenuation measure z/(dr) = n{r) dr, 
n G £j'q(,(M+), then the dispersion- attenuation density n can be calculated by 
applying the formula for inversion of the Stieltjes transform: 

n{r) = — lim Iing{—r — ie) (20) 

7t £-s>0+ 

|42j . More generally, Imf7(p) > if Imp < 0. The distributional limit 
lim£^o+ Im5f(— r — ie) is non-negative, hence it is a non-negative measure. 
For a general non-negative measure z/ we thus have the inversion formula 

i/([a, 6])= / Im lim g{—r — ie) (21) 

for an arbitrary interval [a, 6] of the real axis. 



Eq. (20) can be used to determine the dispersion-attenuation measure 



when the dispersion-attenuation function is given in the form of an analytic 



function. On the other hand either eq. (28) or eq. (29) can be used to deter- 
mine the spectral density hi. This function is subject to only two constraints: 
(1) it is non-negative; (2) it must decrease faster that l/ln(r) as r — )■ oo. 

The Green's function can be calculated by an inverse Laplace transform 
which involves integration over the imaginary axis p = —iu. It is therefore 
important to investigate the functions b, V and A on the imaginary axis in 
the p plane. Substituting p = — iw in the integral representation of b we get 
the following frequency-domain formulas: 

Aiiu) := Rebi-iu) = u^ f ^^ (22) 

Vi{uj) := -\mb{-iuj) =uj f ^^|^ (23) 

i]0,oo[ W + ^ 

For completeness we also note that the propagation speed c{u) and the 
quahty factor Q{uj) can be expressed in terms of the measure v 

^ = -Imi?(-ia;)/a; ^ 1 + M^ (24) 

c[uj) Co OJ 
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It can now be proved that the attenuation is subhnear as a function of 
real frequency. 

Theorem 3.1 

As a; —7- oo 

1. Ai{ijj) = o[uj] and Vii^uj) = o[u]; 

2. ifM := jjp^j z/(dr) < cx) then Ai{uj) -^ M. 
Proof. For cu > 1 

1 u u 1 u 1 

2 u;2 + r^ (w + r)2 l + rw + r 1 + r 

The integrand of Ai{u))/uj is thus bounded by an integrable function and it 
tends to for w — t- oo. The thesis follows from the Lebesgue Dominated 
Convergence Theorem and eq. (tsl). 

For M = f.Q, z/(dr) < oo we note that the integrand of ^1(0;) tends to 
1 and thus ^1(0;) — )■ M by the Lebesgue Dominated Convergence Theorem. 

For oj > 1 the integrand r/i^u"^ + r^) of Vi{u)/u is bounded by an inte- 
grable function 

1 r r r 1 1 

< -. — < < 



2 cj^ + r^ (u + ry 1 + rl + r 1 + r 
Since the integrand tends to for w — )■ 00, Vi{uj)/uj = o[l]. 

By similar arguments one can prove that for u; — > 



n 



Ai{uj)=o[u^] (26) 

Vi{uj) = o [io^] (27) 



Since x/{l + x) is an increasing function, eq. (22) implies that the atten- 
uation function ^1(0;) is non-decreasing function of u; > 0. 

If z/ is absolutely continuous with respect to the logarithmic scale on the 
r axis then there is a locally integrable function h{r) such that z/(dr) = 

11 



h{r)dr/r and L [h{r)/{l + r)]dr/r < oo. The function b now assumes the 
form b{p) = P Jq h{r)/[p + r] dr. The variable r can be replaced by a variable 
r = 1/r with the dimension of time. Let hi{T) := /i(l/r). Substituting 
p = —iu as well we get a formula used in physical acoustics 



2 / ^hiir) 



and its companion formula 

where hi > and 

-^^^dr < cx) 30 

1 + r ^ ^ 

The variable r is different from the relaxation time r = 1/r appearing in 
eq. ([3|. For example in the case of a finite relaxation spectrum r„ = l/r„ 
(eq. ^5|) the attenuation spectrum is also discrete but different from the 
relaxation spectrum. The attenuation spectrum can be calculated from the 
relaxation spectrum by solving a set of non-linear equations following from 
the equation 

Qip)B{pf=pp' (31) 

and conversely. 

4. Implications of finite speed of wave propagation. 

Since 1/Q{p) = pJ{p) = J(0) + J'{p), where J is the creep comphance, 



the Green's function (15) can be expressed in the form of a convolution 



g{t,x) = J'{t)*H{t- \x\/co,\x\) + JoH{t- |x|/co, |x|) 

where H{T,r) is the inverse Laplace transform of the function e~^^P'*^'/(47tr) 
In terms of the inverse Fourier transformation 



^7r2 r ., _^ 



oo 



H(T,r) = ^^ e~'^"-''(-''^)"dw 



Note that b(—iu) = b{iu). The integrand is square integrable if r > and 

/■oo 

/ e-2^*'''(-^^)"dw<oo (32) 

Jo 
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0.005 



Figure 1: Two snapshots of Green's functions for power law attenuation with a = 3/2 
(soHd Une) and a = 1/2 (dashed hne). 



If eq. (32) holds then the Paley- Wiener theorem (Theorem XII in j43j ) 
can be apphed. By this theorem H{t, r) = for r < if and only if 



Re6(— icj) 
The function J' is causal. Hence, if H{t, r 

ft-lxl/co 



dw < oo 

for r < then 



(33) 



vanishes for t < \x\/co. 



J'{s) H{t — \x\/co — s, \x\) ds 







In particular, if b{p) ~oo ap", a > 0, Rep > 0, the inequalities (32) and 



(33) are satisfied if < a < 1. 



It follows from Theorem 3.1 that the asymptotic growth of of the atten- 
uation function Re6(— icj) in the high-frequency range is sublinear. This is 



however insufficient to satisfy inequality (33) and vanishing of the wave field 
ahead of the wavefront |x| = cot. For example for b{—iu)) ~oo w/ln"(cj) 



the integral in (33) does not always converge. Since we are interested in the 
high-frequency behavior, we can replace 1 + u"^ hj w^ in the denominator of 
(p|. For a = 1 

dy 



y My) 



ln(ln(w)) 
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is unbounded for cj — )• oo, hence (33) is not satisfied. For a 7^ 1 



dy ln^-"(w) 



y ln"(y) 1 - a 

is unbounded for a < 1 and bounded for a > 1. We thus see that in contrast 
to the power law attenuation, hnear growth is not a hmit case for finite 
propagation speed. Some subhnear cases will also exhibit precursors ahead 
of the wave front. 

In particular the function b{p) = p/ln"(l + p) with < a < 1 is a CBF 
with the asymptotic properties discussed in the previous paragraph. Indeed, 
for Imp > the argument ip = arg(l +p) satisfies the inequality < ip < n. 
Hence ln(l +p) maps the upper half plane into itself and is non-negative for 



p > 0. Hence, by Theorem 2.1, ln(l + p) is a CBF. The same is true for 



ln°(l + p) if < a < 1. In view of a property of CBF functions mentioned 



after Theorem 2.1 this implies that b{p) is a CBF if < a < 1. Moreover 
&(0) = and \im.p^oob{p)/P = O5 hence b{p) is an admissible dispersion- 
attenuation function. We thus have produced an example of an admissible 
dispersion-attenuation function b{p) such that 6(— iw) has sublinear growth 



in the high frequency range but (33) is not satisfied. 



5. Examples of dispersion-attenuation functions. 

5.1. Power-law attenuation. 

Power-law attenuation is commonly used to match experimental disper- 
sion and attenuation data for a wide variety of real viscoelastic materials 
such as polymers, bio-tissues and some viscoelastic fluids. 

Consider the viscoelastic medium defined by the following equation of 
motion 

p {Dl + 2a DJ.+" + a^ D^") u = AV^ u + 5{t) 6{x) (34) 

with the initial data u{0,x) = M,f(0,x) = in d dimensions, d = 1,3. D^ 
denotes the Caputo fractional derivative of order a 

DS/(^) := r fe^/'(^) d^ (35) 

Jo Ul - a) 
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It is assumed that y4>0, a>0, 0<a<l. The Laplace transformation 
Ct with respect to the time variable and the Fourier transformation J-"^ with 



respect to the spatial variable bring eq. (34) to the following form 
pgipY u{p, k) = —A k^ u{p, k) + 1 



(36) 



where u := T^ {Ct{u)) and g{p) := p + ap"". Eq. (36) implies that in the case 
of power law attenuation Q{p) = pA/g{py and B{p) = g{p)/co. 

We shall begin with solving eq. (p4J) in one dimension. Applying the 



inverse Fourier transformation to eq. (|36|) 

u^^\p, x) 



-^ikx+pt 



-dk 



271 J_^pg{py + Ak 

The contour can be closed by a large half-circle in upper-half complex A;-plane 
if X > 0, in the lower- half complex fc-plane if x < and 



1 



Ak^ + pg{p)^ 2ig{p)A 



1 



k-ig{p)/co k + ig{p)/co 



where Cq := ^ Aj p. We now restrict ourselves to imaginary values of p = 
—\w. Since Im[i5f(— iw)] = ^Qg{—\vo) = a\w\°' cos((l — a)n/2) > 0, the 
residuum at ±ig{p)/co contributes if ±x > 0. Hence 



u{p,x) 



I Q-9(.P)\x\/co+pt 



2A 



g{p) 



and 



M«(t,x) 



1 



^pt-g{p)\x\/co 



dp 



(37) 



47iiAj_i^g{p) 

The solution u^^'^ of (34) in a three-dimensional space is given by the formula 

-1 d 



u^^\t,x) 



2nr dr 



u^'\t,r) 



Hence 

u^^\t,x) 



1 



87r2i|x| CqA 



_,p(i-|x|/co)-ap" \x\ 



dp 



4na^/° |x|^+i/"cov4 



Pa{{t-\x\/co)/{a\x\r) (38) 
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1 ^ico 

Po.{z):=— ey'e-y"dy (39) 

27n J_i^ 



where 



The function P^ is a totally skewed Levy stable probability density 

For a = 1/2 the a-stable probability Pa can be expressed in terms of 
elementary functions. For a = 1/3, 2/3 the a-stable probability can be ex- 
pressed in terms of Airy functions, see Ref. [17]. The Green's functions are 
easy to calculate explicitly in these cases. The Green's function vanishes at 
and outside the wavefront, as can be seen in the case al/2 in Fig. fl] For 
a > 1 there is no wavefront and the signal peak is preceded by a precursor 
extending to infinity, as can be seen in the case of a = 3/2 in Fig. [1} The 
limit case, as we shall see, is the asymptotic behavior b{p) ~oo ap/\n{p). In 
the limit case the propagation speed is unbounded. 

Since Re6(— iw) = aw° cos(7ra/2) > and for p = —iu 



|gPi-g(p)r/co| _ Q-o-'^" cos(a7t/2) 



the integrals 

1 r°° 

^^^'^^ ^^ 2^y (-^TP''p"'9ipr'''e^'''^P^^/''^dp 

are uniformly convergent for r > and a theorem in [18], Sec. 4.44 implies 
that the derivatives (9"+™- f/dt"' dr"^ exist for all positive integers n, m and 
r > 0. For < a < 1 the Green's function vanishes outside the wavefront and 
is infinitely different iable everywhere. Hence the Green's function vanishes 
at the wavefront with all its derivatives |7F] . 

We shall conclude the section with a spectral characterization of the power 
law attenuation: 

Theorem 5.1 

The spectral measure of the dispersion- attenuation function ap", < a < 1, 
is 

„(do = <.!i^^<^r-'d« (40) 

n 
Proof. Substituting the Laplace transform 

y-" = / e-y' [z^-^/Tia]] dz (41) 

Jo 
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in the Laplace transform 

x°-i = / e-'^y [i/-°/r(l - a)] dy (42) 

yields the Stieltjes transform 



X = / az 



n Jq X + z 



which implies eq. (40). 



n 



5.2. Eq. (14) does not imply that the relaxation modulus is completely mono- 
tone. 

All the properties of the dispersion-attenuation function b derived in Sec. [3] 
follow from the fact that B{p) = p p^^'^ / Q^pY^'^ is a complete Bernstein func- 
tion. It follows from the above property of B(p) that Q(p)^/^ is a complete 
Bernstein function [H]. On the other hand Q need not be a complete Bern- 
stein function. Therefore G need not be a completely monotone function 
even though i? is a complete Bernstein function. 

An example is provided by the power law attenuation model. The func- 
tion i? is a complete Bernstein function for < a < 1. On the other hand Q 
is a complete Bernstein function and G is completely monotone if and only 



if 1/ < a < 1. (Fig. ^. Another example is given in Section 5.4 



5.3. Cytoskeleton rheology. 

The continuum model of cytoskeleton rheology [IHl EI] is based on a 
simple creep compliance in pure extension J(t) = Jit^/T{l + (3) with Ji > 
0, < /3 < 1. The last inequality follows from the fact that the creep 
compliance in a viscoelastic material with a positive relaxation time spectrum 
is a Bernstein function[5j. A typical value is /3 = 0.24. It follows that 
Q[p) = Jf V, Co = oo, B{p) = b{p) = {pJif/^p'^ with a := 1 - /3/2 < 1. 
This model is also confirmed by oscillatory experiments which closely match 
the storage modulus G'{uj) = J{^ co^ cos{n(3/2) and loss modulus G"{uj) = 
J^^u^ sin(7r/3/2) at w ^ 1 Hz. In this case G{t) = Git-f^/T{l - l3), where 
Gi = l/Ji. 



17 



0.5 
0.4 
0.3 
0.2 
0.1 

-0.1 




V 



V/ 



20 



40 



60 



80 



100 



Figure 2: The relaxation modulus G for the power law attenuation. The exponent values 
are a = 0.3, 0.4, 0.5, 0.6, 0.9, from bottom to top. 



The fractional Kelvin- Voigt model of bio-tissue 
relaxation modulus 

G{t) = G^ + G^t-P/T{l-P) 



is defined by the 



(43) 



with Goo, Gi > and < /3 < 1. The last inequality follows from the fact 
that the relaxation modulus of a viscoelastic medium with a positive is a 
completely monotone function. In this case Q{j)) = Goo + Gip^, Cq = oo 
and B{p) = b{p) ~o [p/Gi] ^''^ p", where a := 1 - /3/2 < 1. For t ^ the 
fractional Kelvin- Voigt relaxation modulus is asymptotically equivalent to 
the relaxation modulus of the Rouse model of dilute polymer solutions [7j 
and to the relaxation modulus of a derived in |51] for a ladder model of soft 
tissue. 

In all these cases the propagation speed is unbounded because of an infi- 
nite value of the relaxation modulus at t = 0. In all these cases the attenua- 
tion function follows a subhnear power law. 



5.J^-. Discrete dispersion-attenuation spectra. 

If the measure z/ = J2n=i (^n£r„, Tn, Cn > 0, where Ea 
the Dirac measure concentrated at the point a, then 

N 

pJ2 



h{p) 



Cn 



n=l 



p 



Tr. 



5{r — a) denotes 



(44) 



It is however possible to prove that (44) imphes that Q is not a complete 



Bernstein function. The proof is not very short although quite elementary. 



Therefore eq. (44) implies that the relaxation modulus G is not completely 
monotone and the relaxation time spectrum is not positive. It follows that 
the dispersion-attenuation spectra of viscoelastic media always have a non- 
trivial continuous component. 

6. Superlinear power law attenuation and the K-K dispersion re- 
lations. 

A large number of papers in ultrasound acoustics deal with superlinear 
power law attenuation [211 EHl [2ni EH E2] • In order to express attenuation 
in terms of measured dispersion they resort to the K-K dispersion relations 
for the wave number K{u) := iB{—iu). It is also assumed that a model 
of acoustic propagation is correct (causal) if the K-K relations with some 
number of subtractions are satisfied. 

The acoustic K-K dispersion relations follow from the assumption that 
the reduced wavenumber function K{u) — oj/cq is the Fourier transform of 
a causal function (or, more generally, causal tempered distribution) L{t). 
A priori the function L{t) has no physical meaning and the assumption of 
causality of L is unwarranted. A physical meaning can be given to the 
function L by the assumption that the equation of motion can be expressed 
in the following form: 

Cq"^ u,tt + (2/co) L*u^t + L*L*u = u^^^ (45) 

Here 

L{p) = B{p) -p/co 

Causality of the integro-differential operator requires that L is a causal func- 



tion. In this case K{uj) = uj/co+iL{—iui) = iuj/co+L{oi}). Eq. (45) is however 
incompatible with the viscoelastic constitutive equations. In a viscoelastic 
equation of motion integral operators should act on the Laplacian of u. This 
kind of causality does not have any physical meaning. It is more important 
that the Green's function is non-causal in the superlinear case. 

In [26] the authors try to guess the viscoelastic constitutive equation 



compatible with ( 45 ) . Their approach involves an approximation of a spatial 



derivative by a temporal derivative. It is however possible to avoid an ap- 



proximation by shifting the dispersive terms on the left-hand side of (45) to 
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the right-hand side. The Laplace transform of the left-hand side of eq. (45) 
is 



P 



l + CoL{p)/p /c, 



-0 



assuming that m(0, x) = u,t(0, x) = 0. Hence (45) has the form 
where G is the inverse Laplace transform of 

r ~ 1 ~^ 

p~^ l + coL{p)/p 



(46) 



Expression (|46|) is the Laplace transform of a completely monotone function 
1 + coL{p)/p /p is the Laplace transform of a Bernstein 



if and only if 

function/ [5lll2|. 

For L{p) = ap'^ 
2cq at^j^"" /T [1 — a) + Cq a?' f^^^ "Yr(3 — 2a). It is obvious that / is a Bernstein 
function (and G is completely monotone) if and only if 1/2 < a < L We also 
note that L{t) = at^°^^/r(— a) is a distribution. Convolution with L(t) is 
a Riemann-Liouville fractional differential operator of order a. The order of 



with a,a > 0, the function / assumes the form 6{t) + 



the fractional differential equation (45) is 2 if a < 1 and 2q; if a > 1. (In the 



literature the highest-order derivative L*L*u is often incorrectly neglected). 
For a > 1 the order of the time-differential operator exceeds the order of the 
spatial differential operator, hence the equation is formally parabolic. 

Superlinear power law attenuation results in a rather strange behavior of 
the dispersion function. Note that the attenuation function 



A{uj) := ReB{—iu) = a cos{an/2) \ijj\ 



(47) 



is non-negative. Hence, if < a < 1 then the constant a must be non- 
negative. On the other hand, for 1 < a < 3 the constant a must be 
non-positive. As the frequency tends to infinity the frequency-dependent 
phase speed c{u), given by the equation l/c(w) := Re[iB{—iijj)]/u = 1/cq + 
a |ci;|"~^ sin(a7r/2) increases to its maximum value Cq if < a < 1 ("ab- 
normal dispersion"). For 1 < a < 2 it increases from cq at zero frequency 
to infinity at a finite frequency coi = [cq |a| sin(a7t/2)] '^°~ ' and changes 
sign. For 2 < a < 3 phase speed decreases from cq at a; = to at infinite 
frequency ("normal dispersion"). This behavior of the phase speed results in 
the appearance of precursors. 
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7. Viscoelastic constitutive equations with a Newtonian viscosity 
component. 

We shall briefly consider viscoelastic constitutive equations with a New- 
tonian component: 

CT(t) = f G{t- s) e{s) ds + r] e{t) (48) 

where the Newtonian viscosity r/ > 0. If G is completely monotone then the 
function Q{p) = pG{p) +rip is a complete Bernstein function and the theory 
developed above remains valid except that limp__i.oo Q{p) = oo. 

A particular case is the hysteretic damping model deflned by the consti- 



tutive equation (48) and eq. (43) 



Since the limits limp_^o [p G{p)] = Goo and linip^oo [pGip)] = Gq are flnite, 
5(0) = and B{p) ~oo (p/vY^^P^^^- Thus c = and B{p) = b{p). This 
implies that the Green's function does not have any wave fronts. Furthermore 
it follows from the theory of asymptotics of Stieltjes transforms that the 
dispersion-attenuation measure has an asymptotic behavior consistent with 
the asymptotics of b{p): 

j,(dr)~o^^^^r-i/Mr (49) 

Hence high frequency asymptotics of the attenuation and dispersion functions 
is entirely determined by the Newtonian component: 

Mco) -oo - (P/V)'^'^' / ^— ^dr = - ip/vY^'io'/' (50) 

7t /n r"= + Co'"' 7r 



Similarly 



^iM -o^lip/vY^'u: /^^dr = A(p/^)i/2^i/2 (51) 

n J r"^ + 00"^ 6 71 



The low frequency asymptotics of the attenuation and dispersion is en- 
tirely different. For small p the term pG{p) ~o Goo > 

1/0 

B{p) -^oP I'^ + VP/Goo] /coo ~oP [1 -'7P/(2Goo)]/coo 



where Coo := a/Gq/p- Hence 

B{—iijj) = —iu/coo + f3{—iuj) 
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where 



/3(-M 



V 



2Gr 



-u 



(52) 



This shows that in the low frequency range the attenuation function has an 
approximately quadratic frequency dependence. 

The above described low- and high-frequency behavior described here has 
been experimentally confirmed for superionic glasses, with the high-frequency 
asymptotic exponent a ~ 0.8 ^20j . 



The inverse Fourier transform (12) involves the exponential 



g~i(j(t-r/coo) g- 



-iui) r 



For large t — r/coo eq. (52) can be substituted in (12). Upon the changing 
the integration variable y = to r^/^ and setting Q{p) ~ Q(0) = Goo in the 



integrand eq. (12) assumes the form 



^(t,x) = /((t-r/Coo)/ri/2)/r3/2 



(53) 



where 



/(^) 



87t2a 



^-iyz^-r,y'/(2G 



OO L-OO 



dy 



V2c^' 



87r3/2 cr V2 ^1/2 



-riz'^/{8GooCoo) 



(54) 
According to eq. (53 ) the Green's function in the long-range assumes the form 



of a Gaussian pulse traveling with the speed Coo . The quadratic dependence 
of the attenuation can be associated with the r^/^ scaling of the pulse width. 
A rigorous theory of long-range asymptotics of viscoelastic Green's func- 
tions in the absence of Newtonian viscosity is presented in [S2j. In these 
cases the the pulse form and the scaling of the pulse width are different. 



8. Experimentally measured superlinear power law attenuation in 
an intermediate frequency range. 

The restrictions imposed on the dissipation-attenuation exponent a by 
the positive relaxation spectrum and by finite propagation speed apply to 
the asymptotic behavior of b{p) at infinity. At low frequencies the exponent 
is close to 2. It is likely that the experimentally observed superlinear power 
law behavior applies to an intermediate frequency range. 
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For limited frequency bands very good fits of viscoelastic models to su- 
perlinear power laws can be obtained. For example, the Voigt model 

a = Gqc + Tje (55) 

is a special case of the class of viscoelastic models considered here. In this 

case 



puj"^ 






c{u) 



2{GS + v'uj^) ^^^^ 



\p(g, + v^gT+F^) 



Both Ai{ijj) and c(a;) have a w^/^ behavior in the high-frequency range, while 
Ai{uj) ~o [p^'^v] I 2 Go^/^j u?, c{uj) ~o Co for w ^ 0. 

In [52] exponential power law attenuation was matched with Voigt models 
in the frequency range 50-600 Hz. For p = 1000 kg m~^. Go = 2.5 and 
lOkPa, rj = 0.5, 2, 4 Pas a reasonable fit was reported with exponent values 
close to 1.5. It is however pointed out in |53] that the frequency range 
covered in data acquisition available with current technology is not adequate 
for demonstrating that the attenuation function satisfies a power law with 
an exponent in the range between 1 and 2. 

Let (3{u) := ln[^i(w)/^i(l)]/ln(u;) denote the local exponent of the at- 
tenuation function. /3 is an analytic function. Consequently if I3{uj) is con- 
stant over a finite range of frequencies then it is constant over the entire fre- 
quency range < a; < oo and the attenuation function Ai satisfies a power 
law. In Figs [3]-[5] the local exponent has been plotted for the data of Ref . [52] 
for the frequency range 50^600 Hz as well as for the corresponding kilo-Hertz 
and megahertz ranges. Approximately constant exponents are observed over 
frequency ranges of logarithmic width ln(50). The solid line represents the 
exponent for Go = lOkPa and rj = 0.5 Pas, the dashed line corresponds to 
Go = lOkPa and ?7 = 4Pas, and Gq = 2.5 kPa and rj = 0.5 Pas. 

9. Conclusions. 

The attenuation and dispersion functions of a viscoelastic medium with 
positive relaxation spectrum have integral representations in terms of a pos- 
itive Radon measure. Each positive Radon measure satisfying the growth 
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Figure 3: Local exponent in the frequency range 4- 600 Hz. 
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Figure 4: Local exponent in the frequency range 1 -^ 6 kHz. 
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Figure 5: Local exponent in the frequency range 1 -;- 6 MHz. 



condition rtSj) determines an admissible attenuation function and the cor- 
responding dispersion function. The attenuation function has a subhnear 
attenuation growth for a; -t- oo. 

Superhnear growth of logarithmic attenuation rate for a; — ?■ oo is in- 
compatible with hereditary viscoelastic constitutive equations with positive 
relaxation spectrum. It is also incompatible with finite propagation speed. 
Superhnear growth of attenuation results in a precursor of infinite extent 
preceding the main signal. 

Superhnear dependence of attenuation on frequency observed in some 
viscous materials is possible in the intermediate frequency range. 
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